# step10 Cluster the cells
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)

table(pbmc$seurat_clusters)








() FindClusters()
#$ find . | grep "R$" | xargs grep -n "^FindClusters" --color=auto 2>/dev/null

#./seurat-4.1.0/R/generics.R:76:FindClusters <- function(object, ...) {

聚类的确定

#' Cluster Determination
#'
#' Identify clusters of cells by a shared nearest neighbor (SNN) modularity
#' optimization based clustering algorithm. # 通过 SNN 的模块优化，
#' First calculate k-nearest neighbors and construct the SNN graph. #1.计算 KNN 和 SNN
#' Then optimize the modularity function to determine clusters. #2.优化模块化方法来确定聚类结果。
#'
#' For a full description of the algorithms, see Waltman and #完整论文见这一篇
#' van Eck (2013) \emph{The European Physical Journal B}. Thanks to Nigel
#' Delaney (evolvedmicrobe@github) for the rewrite of the Java modularity 
#' optimizer code in Rcpp!
#' ND这个人又重写了一个Java版的，可能是: https://github.com/evolvedmicrobe/seurat/tree/master/java
#'
#' To run Leiden algorithm, you must first install the leidenalg python
#' package (e.g. via pip install leidenalg), see Traag et al (2018). 
#' #如果运行 Leiden 算法，需要依赖pyhton包 leidenalg
#'
#' @param object An object
#' @param ... Arguments passed to other methods
#'
#' @return Returns a Seurat object where the idents have been updated with new cluster info;
#' latest clustering results will be stored in object metadata under 'seurat_clusters'.
#' Note that 'seurat_clusters' will be overwritten everytime FindClusters is run
#' # 返回值是 Seurat 对象，每次运行本函数都覆盖更新 'seurat_clusters'
#'
#' @export
#'
#' @rdname FindClusters
#' @export FindClusters
#'
FindClusters <- function(object, ...) {
  UseMethod(generic = 'FindClusters', object = object)
}








()FindClusters.Seurat
#./seurat-4.1.0/R/clustering.R:414:FindClusters.Seurat <- function(

#' @importFrom methods is
#'
#' @param graph.name Name of graph to use for the clustering algorithm
#'
#' @rdname FindClusters
#' @export
#' @concept clustering
#' @method FindClusters Seurat
#'
FindClusters.Seurat <- function(
  object,
  graph.name = NULL,
  modularity.fxn = 1, #该参数干啥的？
  initial.membership = NULL,
  node.sizes = NULL,
  resolution = 0.8,  #分辨率，重要参数。默认 0.8。我们传入自定义的0.5
  method = "matrix",
  algorithm = 1, #默认算法1
  n.start = 10,
  n.iter = 10,
  random.seed = 0, #随机数种子
  group.singletons = TRUE, #分组的奇异值怎么处理？ 如果是F，则分到一组名叫 sigletion
  temp.file.location = NULL, #临时文件位置
  edge.file.name = NULL, #边文件名
  verbose = TRUE, #输出过程
  ...
) {
  # (A1) 参数检查
  CheckDots(...)

  # (A2) 传入为空，则使用默认值 RNA_snn
  graph.name <- graph.name %||% paste0(DefaultAssay(object = object), "_snn")

  # 检查该值，如果不在 对象中，则报错
  # > names(pbmc)
  #[1] "RNA"     "RNA_nn"  "RNA_snn" "pca"     "umap"    "tsne"   
  if (!graph.name %in% names(x = object)) {
    stop("Provided graph.name not present in Seurat object")
  }
  # 对象中的 object[[graph.name]] 是不是 Graph 类？如果不是，则报错
  if (!is(object = object[[graph.name]], class2 = "Graph")) {
    stop("Provided graph.name does not correspond to a graph object.")
  }

  # (A3) 甩给 FindClusters.default 处理，返回一列df，带列名
  clustering.results <- FindClusters(
    object = object[[graph.name]], #这里是一个 Graph 对象，其实是稀疏矩阵的子类，就是 RNA_snn 矩阵
    modularity.fxn = modularity.fxn,
    initial.membership = initial.membership, #default null
    node.sizes = node.sizes, #default null
    resolution = resolution,
    method = method,
    algorithm = algorithm, #默认1
    n.start = n.start,
    n.iter = n.iter,
    random.seed = random.seed, #随机数种子 0
    group.singletons = group.singletons, #默认 T
    temp.file.location = temp.file.location, #默认 null
    edge.file.name = edge.file.name, #默认 null
    verbose = verbose, #默认T
    ...
  )

  # 拼凑列名: "RNA_snn_res.0.5"
  colnames(x = clustering.results) <- paste0(graph.name, "_", colnames(x = clustering.results))

  #(A4) 添加到 meta.data
  object <- AddMetaData(object = object, metadata = clustering.results)
  
  #(A5) 更新 Idents，列名(其值只有一列)，取第一列
  Idents(object = object) <- colnames(x = clustering.results)[ncol(x = clustering.results)]

  #(A6)提取水平
  levels <- levels(x = object)

  #(A7) 把水平强制变为数字
  # 遇到 warning 和 error 则直接返回原 水平
  levels <- tryCatch(
    expr = as.numeric(x = levels),
    warning = function(...) {
      return(levels)
    },
    error = function(...) {
      return(levels)
    }
  )

  #(A8) 给 Idents 重新赋值：因子，及其水平
  Idents(object = object) <- factor(x = Idents(object = object), levels = sort(x = levels))

  #(A9) 更新 meta.data$seurat_clusters 为 Idents
  object[['seurat_clusters']] <- Idents(object = object)



  #(A10) 记录日志
  cmd <- LogSeuratCommand(object = object, return.command = TRUE)
  #更新 cmd@assay.used
  slot(object = cmd, name = 'assay.used') <- DefaultAssay(object = object[[graph.name]])
  # 记录到 obj
  object[[slot(object = cmd, name = 'name')]] <- cmd

  # 返回对象
  return(object)
}









() 干活的：根据SNN矩阵，对细胞聚类
#./seurat-4.1.0/R/clustering.R:302:FindClusters.default <- function(

#' @importFrom pbapply pblapply # 带 progress bar 的apply, 速度更慢了！
#' @importFrom future.apply future_lapply #多线程的 lapply
#' @importFrom future nbrOfWorkers # 可用核心个数
#'
#' @param modularity.fxn Modularity function (1 = standard; 2 = alternative).
#' @param initial.membership,node.sizes Parameters to pass to the Python leidenalg function.
#' @param resolution Value of the resolution parameter, use a value above
#' (below) 1.0 if you want to obtain a larger (smaller) number of communities.
#' @param algorithm Algorithm for modularity optimization (
# 1 = original Louvain algorithm; #默认
# 2 = Louvain algorithm with multilevel refinement; 
# 3 = SLM algorithm; 
# 4 = Leiden algorithm). Leiden requires the leidenalg python.
#' @param method Method for running leiden (defaults to matrix which is fast for small datasets).
#' Enable method = "igraph" to avoid casting large data to a dense matrix.
#' @param n.start Number of random starts. 随机起点的数量
#' @param n.iter Maximal number of iterations per random start. 每次随机起点的最大迭代次数。
#' @param random.seed Seed of the random number generator.
#' @param group.singletons Group singletons into nearest cluster. If FALSE, assign all singletons to
#' a "singleton" group #默认T； 如果是F，则把所有奇异值分给一个叫sigletons 放到 ‘singleton’ 中。
#' @param temp.file.location Directory where intermediate files will be written.
#' Specify the ABSOLUTE path. 中间文件 文件夹
#' @param edge.file.name Edge file to use as input for modularity optimizer jar. 
#     边文件，作为模块化优化器jar的输入
#' @param verbose Print output
#'
#' @rdname FindClusters
#' @concept clustering
#' @export
#'
FindClusters.default <- function(
  object, #这里是矩阵了
  modularity.fxn = 1,
  initial.membership = NULL,
  node.sizes = NULL,
  resolution = 0.8,
  method = "matrix",
  algorithm = 1,
  n.start = 10, #随机开始的次数
  n.iter = 10, #每次最大迭代次数
  random.seed = 0,
  group.singletons = TRUE,
  temp.file.location = NULL,
  edge.file.name = NULL,
  verbose = TRUE,
  ...
) {
  #(A1) 检查参数
  CheckDots(...)

  #(A2) 参数检查
  if (is.null(x = object)) { #如果obj为空，则报错。
    stop("Please provide an SNN graph")
  }
  # 转为小写，algorithm 如果是 louvain，就是方法1
  if (tolower(x = algorithm) == "louvain") {
    algorithm <- 1
  }
  # 如果是 leiden 就是方法4
  if (tolower(x = algorithm) == "leiden") {
    algorithm <- 4
  }


  #(A3) 如果线程数超过1，默认不走这里
  if (nbrOfWorkers() > 1) {
    # 
    clustering.results <- future_lapply(
      X = resolution, #遍历分辨率，看来这里支持输入数组
      FUN = function(r) {
        if (algorithm %in% c(1:3)) { #默认算法是1
          ids <- RunModularityClustering(
            SNN = object,
            modularity = modularity.fxn,
            resolution = r,
            algorithm = algorithm,
            n.start = n.start,
            n.iter = n.iter,
            random.seed = random.seed,
            print.output = verbose,
            temp.file.location = temp.file.location,
            edge.file.name = edge.file.name
          )
        } else if (algorithm == 4) {
          ids <- RunLeiden(
            object = object,
            method = method,
            partition.type = "RBConfigurationVertexPartition",
            initial.membership = initial.membership,
            node.sizes = node.sizes,
            resolution.parameter = r,
            random.seed = random.seed,
            n.iter = n.iter
          )
        } else {
          stop("algorithm not recognised, please specify as an integer or string")
        }
        # 分组结果加名字，就是原 cid
        names(x = ids) <- colnames(x = object)
        
        # 聚类奇异值?
        ids <- GroupSingletons(ids = ids, SNN = object, verbose = verbose)
        
        # 结果变因子
        results <- list(factor(x = ids))
        # 加列名， res.+分辨率
        names(x = results) <- paste0('res.', r)
        # 返回结果
        return(results)
      }
    )

    # 强制转为数据框
    clustering.results <- as.data.frame(x = clustering.results)
  

  #(A3) 如果是单线程 (默认，是不开多线程的)
  } else {
  	# 创建空字符串，只有行名，其他为空。
    clustering.results <- data.frame(row.names = colnames(x = object))

    # 遍历 分辨率
    for (r in resolution) {
      # 默认算法是1
      if (algorithm %in% c(1:3)) {
        ids <- RunModularityClustering(
          SNN = object,
          modularity = modularity.fxn,
          resolution = r,
          algorithm = algorithm,
          n.start = n.start,
          n.iter = n.iter,
          random.seed = random.seed,
          print.output = verbose,
          temp.file.location = temp.file.location,
          edge.file.name = edge.file.name)
      } else if (algorithm == 4) {
        ids <- RunLeiden(
          object = object,
          method = method,
          partition.type = "RBConfigurationVertexPartition",
          initial.membership = initial.membership,
          node.sizes = node.sizes,
          resolution.parameter = r,
          random.seed = random.seed,
          n.iter = n.iter
        )
      } else {
        stop("algorithm not recognised, please specify as an integer or string")
      }

      # named vector， 添加的是 cid
      names(x = ids) <- colnames(x = object)

      # 奇异值是否聚类？
      ids <- GroupSingletons(ids = ids, SNN = object, group.singletons = group.singletons, verbose = verbose)

      # 添加一列
      clustering.results[, paste0("res.", r)] <- factor(x = ids)
    }
  }
  return(clustering.results)
}










#



() (java 程序的 C++端口)
$ find . | grep "R$" | xargs grep -n "RunModularityClustering" --color=auto 2>/dev/null
看样子又到头了。R代码的尽头通常是C++，这里又出来了 java。
#./seurat-4.1.0/R/clustering.R:1737:RunModularityClustering <- function( #


# Runs the modularity optimizer (C++ port of java program ModularityOptimizer.jar)
#
# @param SNN SNN matrix to use as input for the clustering algorithms
# @param modularity Modularity function to use in clustering (1 = standard; 2 = alternative)
# @param resolution Value of the resolution parameter, use a value above (below) 1.0 
#    if you want to obtain a larger (smaller) number of communities
# @param algorithm Algorithm for modularity optimization 
# (1 = original Louvain algorithm; 
# 2 = Louvain algorithm with multilevel refinement; 
# 3 = SLM algorithm; 
# 4 = Leiden algorithm). Leiden requires the leidenalg python module.
# @param n.start Number of random starts
# @param n.iter Maximal number of iterations per random start
# @param random.seed Seed of the random number generator
# @param print.output Whether or not to print output to the console
# @param temp.file.location Deprecated and no longer used
# @param edge.file.name Path to edge file to use
#
# @return Seurat object with identities set to the results of the clustering procedure
#
#' @importFrom utils read.table write.table
#
RunModularityClustering <- function(
  SNN = matrix(),
  modularity = 1,
  resolution = 0.8,
  algorithm = 1,
  n.start = 10,
  n.iter = 10,
  random.seed = 0,
  print.output = TRUE,
  temp.file.location = NULL,
  edge.file.name = NULL
) {
  # 不提供则默认为 ""
  edge_file <- edge.file.name %||% ''

  # 调用C++
  clusters <- RunModularityClusteringCpp(
    SNN,
    modularity,
    resolution,
    algorithm,
    n.start,
    n.iter,
    random.seed,
    print.output,
    edge_file
  )
  # 返回聚类结果
  return(clusters)
}







==> 下一节要看看 Rcpp /C++了。被迫中断好几次了。
#./seurat-4.1.0/R/RcppExports.R:4:RunModularityClusteringCpp <- function(SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename) {
#./seurat-4.1.0/R/RcppExports.R:5:    .Call('_Seurat_RunModularityClusteringCpp', PACKAGE = 'Seurat', SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename)

RunModularityClusteringCpp <- function(SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename) {
    .Call('_Seurat_RunModularityClusteringCpp', PACKAGE = 'Seurat', SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename)
}

R代码到 .Call() 这里就是尽头了，接着追踪就是Rcpp了。











() 看 C++文件 src/RcppExports.cpp
这是一个中间文件。可以跳过。
$ find . | grep -v "R$" | xargs grep -n "RunModularityClusteringCpp" --color=auto 2>/dev/null
./seurat-4.1.0/src/RcppExports.cpp:14:// RunModularityClusteringCpp
./seurat-4.1.0/src/RcppExports.cpp:15:IntegerVector RunModularityClusteringCpp(Eigen::SparseMatrix<double> SNN, int modularityFunction, double resolution, int algorithm, int nRandomStarts, int nIterations, int randomSeed, bool printOutput, std::string edgefilename);
./seurat-4.1.0/src/RcppExports.cpp:16:RcppExport SEXP _Seurat_RunModularityClusteringCpp(SEXP SNNSEXP, SEXP modularityFunctionSEXP, SEXP resolutionSEXP, SEXP algorithmSEXP, SEXP nRandomStartsSEXP, SEXP nIterationsSEXP, SEXP randomSeedSEXP, SEXP printOutputSEXP, SEXP edgefilenameSEXP) {
./seurat-4.1.0/src/RcppExports.cpp:29:    rcpp_result_gen = Rcpp::wrap(RunModularityClusteringCpp(SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename));



// RunModularityClusteringCpp
IntegerVector RunModularityClusteringCpp(Eigen::SparseMatrix<double> SNN, int modularityFunction, double resolution, int algorithm, int nRandomStarts, int nIterations, int randomSeed, bool printOutput, std::string edgefilename);
RcppExport SEXP _Seurat_RunModularityClusteringCpp(SEXP SNNSEXP, SEXP modularityFunctionSEXP, SEXP resolutionSEXP, SEXP algorithmSEXP, SEXP nRandomStartsSEXP, SEXP nIterationsSEXP, SEXP randomSeedSEXP, SEXP printOutputSEXP, SEXP edgefilenameSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< Eigen::SparseMatrix<double> >::type SNN(SNNSEXP);
    Rcpp::traits::input_parameter< int >::type modularityFunction(modularityFunctionSEXP);
    Rcpp::traits::input_parameter< double >::type resolution(resolutionSEXP);
    Rcpp::traits::input_parameter< int >::type algorithm(algorithmSEXP);
    Rcpp::traits::input_parameter< int >::type nRandomStarts(nRandomStartsSEXP);
    Rcpp::traits::input_parameter< int >::type nIterations(nIterationsSEXP);
    Rcpp::traits::input_parameter< int >::type randomSeed(randomSeedSEXP);
    Rcpp::traits::input_parameter< bool >::type printOutput(printOutputSEXP);
    Rcpp::traits::input_parameter< std::string >::type edgefilename(edgefilenameSEXP);
    rcpp_result_gen = Rcpp::wrap(RunModularityClusteringCpp(SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename));
    return rcpp_result_gen;
END_RCPP
}




./seurat-4.1.0/src/RcppExports.cpp:408:    {"_Seurat_RunModularityClusteringCpp", (DL_FUNC) &_Seurat_RunModularityClusteringCpp, 9},

static const R_CallMethodDef CallEntries[] = {
    {"_Seurat_RunModularityClusteringCpp", (DL_FUNC) &_Seurat_RunModularityClusteringCpp, 9},
...





() 看 C++文件 src/RModularityOptimizer.cpp
干活的。输入 SNN 矩阵，返回整数向量（每个细胞的类别数字）。
警告：本函数没有执行和测试，所以其中的注释全靠目测和瞎猜，仅供参考。
我对这个算法还是很感兴趣的，无奈看不懂 C++ 代码。 以后 C++ 过关了接着看。



./seurat-4.1.0/src/RModularityOptimizer.cpp:23:IntegerVector RunModularityClusteringCpp(Eigen::SparseMatrix<double> SNN,


// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppProgress)]]
// [[Rcpp::export]]
IntegerVector RunModularityClusteringCpp(Eigen::SparseMatrix<double> SNN,
    int modularityFunction,
    double resolution,
    int algorithm,
    int nRandomStarts,
    int nIterations,
    int randomSeed,
    bool printOutput,
    std::string edgefilename) {

  // validate arguments 验证参数合法性
  if(modularityFunction != 1 && modularityFunction != 2) //共2个，标准 和 另一个
    stop("Modularity parameter must be equal to 1 or 2.");
  if(algorithm != 1 && algorithm !=2 && algorithm !=3 && algorithm !=4)
    stop("Algorithm for modularity optimization must be 1, 2, 3, or 4");
  if(nRandomStarts < 1)
    stop("Have to have at least one start");
  if(nIterations < 1)
    stop("Need at least one interation");
  if (modularityFunction == 2 && resolution > 1.0)
    stop("error: resolution<1 for alternative modularity");


  try { //这个括号管的特别长，一直到本函数倒数5-6行。
  bool update;
  double modularity, maxModularity, resolution2;
  int i, j;


  //C++代码的输出
  std::string msg = "Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck";
  // 默认要打印
  if (printOutput)
    Rcout << msg << std::endl << std::endl;



  // Load netwrok
  std::shared_ptr<Network> network;
  if(edgefilename != "") { //默认不走这里
    if (printOutput)
      Rcout << "Reading input file..." << std::endl << std::endl;
    try{
      network = readInputFile(edgefilename, modularityFunction);
    } catch(...) {
      stop("Could not parse edge file.");
    }
  } else {  //默认 edgefilename 是空的，默认走这里
    // Load lower triangle 载入下三角矩阵
    int network_size = (SNN.nonZeros() / 2) + 3; //网络尺寸，非零元素个数
    IVector node1;
    IVector node2;
    DVector edgeweights;
    node1.reserve(network_size);
    node2.reserve(network_size);
    edgeweights.reserve(network_size);

    // 遍历， 列>=行 时跳过，那就是只填充下三角矩阵。
    for (int k=0; k < SNN.outerSize(); ++k){
      for (Eigen::SparseMatrix<double>::InnerIterator it(SNN, k); it; ++it){
        if(it.col() >= it.row()){
          continue;
        }
        node1.emplace_back(it.col());
        node2.emplace_back(it.row());
        edgeweights.emplace_back(it.value());
      }
    }

    // 如果全是0，则报错。
    if(node1.size() == 0) {
      stop("Matrix contained no network data.  Check format.");
    }

    //节点个数，是行数和列数中的最大值。
    int nNodes = std::max(SNN.cols(), SNN.rows());
    // 矩阵转为网络
    network = matrixToNetwork(node1, node2, edgeweights, modularityFunction, nNodes);
    Rcpp::checkUserInterrupt(); //检查用户是否打断？
  }


  // 默认要打印
  if (printOutput)
  {
    Rprintf("Number of nodes: %d\n", network->getNNodes()); //节点数
    Rprintf("Number of edges: %d\n", network->getNEdges()); //边数
    Rcout << std::endl;
    Rcout << "Running " <<  ((algorithm == 1) ? "Louvain algorithm" : ((algorithm == 2) ? "Louvain algorithm with multilevel refinement" : "smart local moving algorithm")) << "...";
    Rcout << std::endl;
  }


  // 重新 计算分辨率
  resolution2 = ((modularityFunction == 1) ? (resolution / (2 * network->getTotalEdgeWeight() + network->getTotalEdgeWeightSelfLinks())) : resolution);

  // 计时器，起点时间
  auto beginTime = duration_cast<milliseconds>(system_clock::now().time_since_epoch());
  std::shared_ptr<Clustering> clustering;
  maxModularity = -std::numeric_limits<double>::infinity();
  JavaRandom random(randomSeed); //第一个java随机数？

  // 开始进度条？
  Progress p(nRandomStarts, printOutput);

  // 对i循环，最大是 参数 最大随机次数
  for (i = 0; i < nRandomStarts; i++)
  {
    //if (printOutput && (nRandomStarts > 1))
    //Rprintf("Random start: %d\n", i + 1);

    VOSClusteringTechnique vosClusteringTechnique(network, resolution2); //这是什么？在 .h 中提到过，但没有更多信息。

    j = 0;
    update = true;
    do //开始迭代
    {
      /*if (printOutput && (nIterations > 1))
        Rprintf("Iteration: %d\n", j + 1);
      */

      if (algorithm == 1) //默认是算法1
        update = vosClusteringTechnique.runLouvainAlgorithm(random); //运行模型
      else if (algorithm == 2)
        update = vosClusteringTechnique.runLouvainAlgorithmWithMultilevelRefinement(random);
      else if (algorithm == 3)
        vosClusteringTechnique.runSmartLocalMovingAlgorithm(random);

      j++; //迭代次数计数器

      // 搜这几个 vosClustering 开头的函数，找到了 java 的文档:
      // https://dwslab.github.io/melt/javadoc_latest/de/uni_mannheim/informatik/dws/melt/matching_jena_matchers/multisource/clustering/VOSClusteringTechnique.html
      // 这里看到了类似的代码，java 代码。
      //https://github.com/dwslab/melt/blob/melt-3.2/matching-jena-matchers/src/main/java/de/uni_mannheim/informatik/dws/melt/matching_jena_matchers/multisource/clustering/ComputeErrDegree.java
      modularity = vosClusteringTechnique.calcQualityFunction(); //计算质量函数

      //if (printOutput && (nIterations > 1))
      //  Rprintf("Modularity: %.4f\n", modularity);
      Rcpp::checkUserInterrupt();
    }
    while ((j < nIterations) && update); //循环条件：j<最大迭代次数，且update不等于0



    if (modularity > maxModularity)
    {
      clustering = vosClusteringTechnique.getClustering(); //获取分类
      maxModularity = modularity;
    }

    /*if (printOutput && (nRandomStarts > 1))
    {
      if (nIterations == 1)
        Rprintf("Modularity: %.4f\n", modularity);
      Rcout << std::endl;
    }*/
    p.increment(); //进度条递增
  }


  // 结束时间
  auto endTime = duration_cast<milliseconds>(system_clock::now().time_since_epoch());
  if(clustering == nullptr) {
    stop("Clustering step failed.");
  }

  // 默认要打印
  if (printOutput)
  {
    if (nRandomStarts == 1) //默认起始次数是10，不走这里
    {
      if (nIterations > 1)
        Rcout << std::endl;
      Rprintf("Modularity: %.4f\n", maxModularity);
    }
    else //默认走这里
      Rprintf("Maximum modularity in %d random starts: %.4f\n", nRandomStarts, maxModularity); //这是使用的R函数？
      Rprintf("Number of communities: %d\n", clustering->getNClusters());
      Rprintf("Elapsed time: %d seconds\n", static_cast<int>((endTime - beginTime).count() / 1000.0));
  }




  // Return results
  clustering->orderClustersByNNodes();
  IntegerVector iv(clustering->cluster.cbegin(), clustering->cluster.cend());
  return iv; //返回聚类结果

  } catch(std::exception &ex) {
    forward_exception_to_r(ex);
  } catch(...) {
    ::Rf_error("c++ exception (unknown reason)");
  }

  return IntegerVector(1);
}












() 看 C++文件 src/ModularityOptimizer.cpp
这个文件和上一个的区别是，文件名首字母少了一个R。
我没有发现调用这个文件的代码，仅仅是名字很像，多看了一眼它的说明部分。


这个代码是翻译的 java包 ModularityOptimizer.jar, 来自 http://www.ludowaltman.nl/slm/
用来聚类，用于早期版本的 Seurat。
评论：这句话意思是说现在的 Seurat 不用这个代码了？
// This code is a translation of the Java package, ModularityOptimizer.jar,
// available from http://www.ludowaltman.nl/slm/ that performs clustering and was
// used by an earlier version of Seurat.


翻译时，接口保留，所以虽然编程语言变了，但结果是一样的。
In translating the code, the interface was maintained, such that the although the programming languages changed, the results are identical. 

所以，重写并没有用更地道的C++写，输出和代码都尽可能和原文java一致。
As a consequence, rather than rewriting the program to be more idiomatic Rcpp/C++,
both the output and the code was kept as close to the original Java as possible. 

比如，为了兼容性， Java 随机数使用 C++ 重写了， std::stable_sort 代替 std::sort 来保证输出和 Java 一模一样。
For example, in order to maintain compatibility, the Java random
number generator was implemented in C++, and std::stable_sort was used in
place of std::sort to match the output results of the Java code exactly.


为了保证C++和Java代码返回值完全一致， C++ 代码加入了原始命令的镜像接口。
// In order to test and verify that the C++ and Java code return the same
// results, the C++ code added also includes a mirror of the original command
// line interface. 

这一版的程序可以这样编译:
This version of the program can be compiled with a command such as:
//
//    clang++ -O3 -std=c++11 -DSTANDALONE -Wall -g ModularityOptimizer.cpp

然后和java同样的调用方式，来确认兼容性，或者识别 输入/输出中的问题。
// And can then be called in an identical fashion to the original Java version to
// verify compatibility or identify any issues in the input/output.




再次查找，确认这个文件哪里使用了？
$ find . | xargs grep -n "ModularityOptimizer" --color=auto 2>/dev/null
./seurat-4.1.0/R/clustering.R:1720:# Runs the modularity optimizer (C++ port of java program ModularityOptimizer.jar)


./seurat-4.1.0/src/ModularityOptimizer.cpp:1:// This code is a translation of the Java package, ModularityOptimizer.jar,
./seurat-4.1.0/src/ModularityOptimizer.cpp:18://    clang++ -O3 -std=c++11 -DSTANDALONE -Wall -g ModularityOptimizer.cpp
./seurat-4.1.0/src/ModularityOptimizer.cpp:23:#include "ModularityOptimizer.h"
./seurat-4.1.0/src/ModularityOptimizer.cpp:31:using namespace ModularityOptimizer;
./seurat-4.1.0/src/ModularityOptimizer.cpp:759:std::shared_ptr<Network> ModularityOptimizer::matrixToNetwork(IVector& node1, IVector& node2, DVector& edgeWeight1, int modularityFunction, int nNodes) {
./seurat-4.1.0/src/ModularityOptimizer.cpp:806:std::shared_ptr<Network> ModularityOptimizer::readInputFile(std::string fname, int modularityFunction) {
./seurat-4.1.0/src/ModularityOptimizer.cpp:840:std::vector<std::string> ModularityOptimizer::split(const std::string& s, char delimiter)
./seurat-4.1.0/src/ModularityOptimizer.h:15:namespace ModularityOptimizer {


./seurat-4.1.0/src/RModularityOptimizer.cpp:13:#include "ModularityOptimizer.h"
./seurat-4.1.0/src/RModularityOptimizer.cpp:15:using namespace ModularityOptimizer;
命名空间一致？ 可能是用到的。



由于C++看不懂，产生好几个悬念：
- 这个 ModularityOptimizer.cpp 在当前版本是否还在用？还在用，至少 vos 开头的部分还在用。
- 这个 ModularityOptimizer.jar 在哪里？用过吗？还是早期也是使用的其 C++ 复刻版？
  * http://www.ludowaltman.nl/slm/
  * https://github.com/CWTSLeiden/networkanalysis
- VOSClusteringTechnique 是从哪里来的？这个类/函数怎么定义的？
  * 搜文件，发现都在 ModularityOptimizer.cpp 中定义的。
- VOS Clustering 又是啥？
  * VOS, which is an abbreviation for visualization of similarities. 
  * https://link.springer.com/chapter/10.1007/978-3-540-70981-7_34 (2007)










() GroupSingletons()
./seurat-4.1.0/R/clustering.R:1314:GroupSingletons <- function(ids, SNN, group.singletons = TRUE, verbose = TRUE) {


把单细胞分组，把单细胞分到他们联系最紧密的组内。

# Group single cells that make up their own cluster in with the cluster they are
# most connected to.
#
# @param ids Named vector of cluster ids
# @param SNN SNN graph used in clustering
# @param group.singletons Group singletons into nearest cluster. If FALSE, assign all singletons to
# a "singleton" group
#
# @return Returns Seurat object with all singletons merged with most connected cluster
#
GroupSingletons <- function(ids, SNN, group.singletons = TRUE, verbose = TRUE) {
  #(A1) 识别 奇点
  # identify singletons
  singletons <- c()
  
  # 找到只出现1次的类
  singletons <- names(x = which(x = table(ids) == 1))
  # 求交集。 为什么还要和原文求交集？ 保证顺序是原文顺序。 测试 ids=c(6,2,2,3,3,3,4,5,5);
  singletons <- intersect(x = unique(x = ids), singletons)
  
  # 如果不分组奇点 
  if (!group.singletons) { #默认不走这里
  	# 则奇点属于一组，名叫 sigleton 且直接返回。
    ids[which(ids %in% singletons)] <- "singleton"
    return(ids)
  }

  #(A2) 计算 奇点 和其他类的联系，归类到联系最紧密的类中
  # calculate connectivity of singletons to other clusters, add singleton
  # to cluster it is most connected to
  # 独特的类名，变字符串
  cluster_names <- as.character(x = unique(x = ids))
  
  # 独特的类名 差集 奇点类
  cluster_names <- setdiff(x = cluster_names, y = singletons)

  #(A3) 连接度，初始化为空向量，长度是 其他类的个数
  connectivity <- vector(mode = "numeric", length = length(x = cluster_names))

  #(A4) 向量命名，就是 非奇点类的名字
  names(x = connectivity) <- cluster_names

  #(A5) 新的类向量 = 老的类向量
  new.ids <- ids

  #(A6) 遍历 奇点类
  for (i in singletons) {
  	# 由于 ids 是 named vactor，比较时是按照值，然后获取其 name，也即是cid
  	# 有这些i是 奇点，所以一次只有一个cid
    i.cells <- names(which(ids == i))

    # 遍历 非奇点类
    for (j in cluster_names) {
      # 获取其细胞的id
      j.cells <- names(which(ids == j))
      # SNN矩阵取子集。 行为奇点一个cid，列为其他类，某个具体类的cid
      subSNN <- SNN[i.cells, j.cells]

      # 随机数种子，怎么会影响 WhichCells() 呢 //todo
      set.seed(1) # to match previous behavior, random seed being set in WhichCells
      
      # 如果是对象，表示是数据框？
      if (is.object(x = subSNN)) {
      	# 只有一行，连接性 = 值为1的个数 / 总列数
        connectivity[j] <- sum(subSNN) / ( nrow(x = subSNN) * ncol(x = subSNN))
      } else { #我认为不走这里。
        connectivity[j] <- mean(x = subSNN)
      }
    }

    # 找连接数最大的，去掉 na
    m <- max(connectivity, na.rm = T)
    # 最大值的序号
    mi <- which(x = connectivity == m, arr.ind = TRUE)
    
    # 找该序号的分类名，如果多个，则抽取1个
    closest_cluster <- sample(x = names(x = connectivity[mi]), 1)
    # 这就是 该奇点 最可能的分类
    ids[i.cells] <- closest_cluster
  }

  # 如果 奇点 长度大于0，则打印：n1 个奇点； n2个最终类。
  if (length(x = singletons) > 0 && verbose) {
    message(paste(
      length(x = singletons),
      "singletons identified.",
      length(x = unique(x = ids)),
      "final clusters."
    ))
  }
  
  return(ids)
}

